---
title: "Beta Diversity: Part II: Gradient Analysis in R"
author: "Michal Kowalewski"
date: "August 2, 2018"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Beta gradient analysis is explictly spatial. We will examine a dataset of marine invertebrate samples from coastal North Carolina. First, read and examine live community data and associated locality info.

We will also need the library 'vegan'.

```{r}
library(vegan)
live <- read.csv('Final_All_Taxa_BetaDiv.csv',head=TRUE)
refs <- read.csv('Refs_local.csv')
```

Now, assess beta diversity gradient as a function of depth difference between localities. Bray-Curtis is a popular measure, but you can use other measures as well.

The obvious way (although perhaps not CPU-efficient) is to use nested for loops. Try those first. 

```{r}
# for loop solution
out1 <- NULL
for (i in 1:nrow(live)) {
 for (j in 1:nrow(live)) {
  if (j > i) {
  dist <- 1 - vegdist(rbind(live[i,], live[j,]))
  D.depth <- abs(refs$Depth[i] - refs$Depth[j])
  habitat <- ifelse(refs$Habitat[i] == refs$Habitat[j], 
                    habitat <- 1, habitat <- 0)
  out1 <- rbind(out1, c(dist=dist, depth=D.depth, hab=habitat,
                      min.n=min(sum(live[i,]), sum(live[j,]))))
  }
 }
}
out3 <- data.frame(out1)
```

You may try to utilize a well tested and efficient function such as vegdist to see if you can make your script faster. 

```{r}
# making use of vegdist function
min.nveg <- cbind(designdist(live, method = "A", terms='minimum'),
                  designdist(live, method = "B", terms='minimum'))
min.nveg2 <- apply(min.nveg, 1, min)
habdist <- vegdist(refs$Habitat, 'euclidean')
hab1 <- ifelse(habdist==0, hab1 <- 1, hab1 <- 0)
out2 <- cbind(dist = 1 - vegdist(live),
              depth = vegdist(refs$Depth, 'euclidean'),
              hab = hab1, min.n=min.nveg2)
out4 <- data.frame(out2)
```

If you created multiple ways to do this exercise, check if those methods agree. Say, your outputs are 'out3' and 'out4'. Check if they are the same.

```{r}
sum(out4!=out3) # check if all values agree
```

The final step is to visualize your results. If you generated multiple solutions. It may be best to write a custom functioon for plotting your data.

```{r}
bg.plot <- function(out) {
bubble <- (out$min.n - min(out$min.n))/(max(out$min.n) - min(out$min.n)) + 1
  b1 <- tapply(out$dist, round(out$depth, -0.5), mean)
  b2 <- tapply(out$dist[out$hab==1], round(out$depth[out$hab==1], -0.5), mean)
  b3 <- tapply(out$dist[out$hab==0], round(out$depth[out$hab==0], -0.5), mean)
plot(out[,2], out[,1], ylim=c(0,1), pch=21, cex=bubble, las=1,
     col=c('tan3', 'skyblue')[as.factor(out$hab)],
     bg=adjustcolor(c('tan3', 'skyblue')[as.factor(out$hab)], 0.3),
     xlab=bquote(Delta~'depth [m]'),
     ylab='Bray-Curtis similarity')
   points(as.numeric(names(b1)), b1, type='l', col='white', lwd=5)
   points(as.numeric(names(b1)), b1, type='l', col='black', lwd=2.5)
   points(as.numeric(names(b2)), b2, type='l', col='white', lwd=5)
   points(as.numeric(names(b2)), b2, type='l', col='tan3', lwd=2.5)
   points(as.numeric(names(b3)), b3, type='l', col='white', lwd=5)
   points(as.numeric(names(b3)), b3, type='l', col='skyblue', lwd=2.5)
  legend('topright', lty=1, col=c('black', 'skyblue', 'tan3'),
       c('all', 'between habitats', 'within-habitat'), cex=0.7)
  legend('topleft', pch=21, pt.cex=c(max(bubble), min(bubble)),
       legend=c(paste('largest n =', max(out$min.n)),
                paste('smallest n =', min(out$min.n))))
}
```

Now, plot your outcomes.

#### Solution 1

```{r}
bg.plot(out3)
```

#### Solution 2
```{r}
bg.plot(out4)
```

#### adding turnover  model (exponential fit)

This function has not been tested(use at your own peril!) 
```{r}
decay<-function(x){
  fitH<- x[which(x[,1]*x[,2]>0),] # get rid of 0 values (not allowed for log operations)
  expH<- lm(log(fitH[,2])~ fitH[,1]) # estimate exponent
  xlH<- seq(0.1,15,0.1) # set a series of points for range of depth difference values (for plotting over a given range of x values)
  ylH<- exp(expH$coefficients[2]*xlH+expH$coefficients[1]) # compute predicted y values (bray-curtis values)
  adjr2 <- summary(expH)$adj.r.squared # compute adjusted r-square
  decay<-cbind(xlH,ylH,expH$coefficients[2])
  return(decay)
}
```

Now plot the plot again but add a decay model line
```{r}
bg.plot(out4)
 points(decay(cbind(out4$depth, out4$dist)), type='l', col='white', lwd=4)
 points(decay(cbind(out4$depth, out4$dist)), type='l', col='red1', lwd=2)
```
And then cite packages...

```{r}
citation('vegan')
```

